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Abstract 

We propose an adaptive biasing algorithm aimed at enhancing the sampling of 
multimodal measures by Langevin dynamics. The underlying idea consists in general- 
izing the standard adaptive biasing force method commonly used in conjunction with 
molecular dynamics to handle in a more effective fashion multidimensional reaction 
coordinates. The proposed approach is anticipated to be particularly useful for re- 
action coordinates, the components of which are weakly coupled, as illuminated in a 
mathematical analysis of the long-time convergence of the algorithm. The strength 
as well as the intrinsic limitation of the method are discussed and illustrated in two 
realistic test cases. 



1 Introduction 

Sampling of multimodal measures is a central problem in many scientific areas, such 
as statistical simulations, in particular molecular dynamics, which constitutes the 
primary focus of the present work. One standard approach to deal with such a situation 
consists in resorting to biasing techniques — e.g. importance sampling methods, in 
order to reduce the multimodal nature of the targeted measure. Under these premises, 
the main difficulty is evidently to devise the correct bias. 

One class of methods proposed in the framework of molecular dynamics and which 
has proven to be useful also for a variety of applications |15| are adaptive biasing 
numerical schemes. The underlying idea here consists in designing adaptively the bias 
such that the new targeted measure be uniform along a priori chosen directions. Only 
these directions have to be chosen, but not the precise analytical expression of the bias. 
Examples of such class of approaches include the Wang-Landau algorithm [53], non- 
equilibrium metadynamics |33[ ITT] and the adaptive biasing force (ABF) method |17| 
which will constitute the main thrust of this contribution. The reader is referred 
to |391 Chapter 5] for a general, mathematically-oriented presentation of adaptive 
methods. 

Let us introduce the ABF method. In what follows, the algorithms will be pre- 
sented in the framework of sampling of configurational space and overdamped Langevin 
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dynamics, but generalization to sampling of phase space and standard Langevin dy- 
namics is straightforward, as can be seen in Section [4] The canonical, Boltzmann- 
Gibbs, measure will be considered here: 

dv(q) = Z- 1 cM-PV(q))dq 1 (1) 

where f3 is proportional to the inverse temperature, q € T>, V : T> — > K is the so-called 
potential energy function, which is assumed to be a smooth function in the following, 
Z = J v exp(— j3V (q)) dq and T> = {q, V(q) < 00} is the configurational space. The 
overdamped Langevin dynamics writes: 



dQ t = -VV(Qt) dt + ^JW^dBt, (2) 

where B t is an n-dimensional standard Brownian motion. Under mild conditions on V, 
this dynamics is ergodic with respect to the measure /1, i.e. trajectory or time averages 
converge to canonical averages. 

For a multimodal measure /i, the sampling obtained with Q t is, however, rather 
poor. Indeed, the typical problem with dynamics ([2| is that the process Q t remains 
trapped for long times in some metastable states. The purpose of an adaptive bi- 
asing force is to enhance sampling by subtracting from V a potential such that the 
aforementioned metastable features are eliminated. This relies on an assumed a priori 
knowledge of those "coordinates" that remain trapped — viz. "slow variables" of the 
dynamics, also called "collective variables" or "reaction coordinates". This method can 
thus be seen as an adaptive importance sampling method. 

In the following, emphasis will be put on the case where two slow variables have 
been identified, namely £1 : T> — > T, and £2 : 2? — > T, where T denotes the one- 
dimensional torus. Throughout the present work, £1 and £2 will be assumed to be 
smooth functions such that |V£i| ^ and |V£2| 7^ 0. 

For simplicity, let us consider for the moment simple reaction coordinates: 

V = T n ,Z 1 (x) = x 1 ,&{x) = x 2 , (3) 

where the components of x are referred to as: x = (x±, X2, x 3 ... n ). The case of general 
£i's will be discussed below, in Sections [3] and |4j In this simple framework, one 
standard ABF-like method is |37j : 

2 

dY t = ( - W + J2 r " V£ Q ) {Y t ) dt + ^/2^dB t , 

for a = 1,2, T?( Xl ,x 2 ) = E(d Xa V(Y t ) | (£i,6)(*t) = (x u x 2 )). 



(4) 



In practice, the dynamics is discretized in time, and the biasing forces, Tf, are approx- 
imated by empirical or time averages in each cells of a grid of the values of (£1 , £2)- The 
bottom line of the method is to observe that, in the long-time limit, (Tj , ) converges 
to VA, where A is the so-called free energy associated to V and (£1,^2) — see [3"5] . 
At equilibrium, the potential eventually becomes V — A o (£i,£a), which is such that 
the associated Boltzmann-Gibbs measure, proportional to exp(— j3(V — Ao (£i,£2))), 
has uniform marginal laws along (£1,^2)- Indeed, the free energy is defined, up to an 
additive constant, as: 

A(x 1 ,x 2 ) = -/3 _1 ln ^ exp(-pv(xi,x 2 ,x 3 ... n ))dx 3 ...n \ . (5) 

Using the definition (|5| of the free energy A, it is easy to check that the equilibrium 
probability density for Q, namely 

tpoo(xi,x 2 ,x 3 ... n ) oc exp ( - /3(V(xi 1 x 2 ,x 3 ... n ) - A(x 1 ,x 2 ))J 

is such that the marginal law of i/'oo along (£i,£2) is uniform: 

1pea(xi,X 2 ,X 3 ... n )dx 3 ... n = 1 T 2 . 

J-n-2 
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Another related important property of the stochastic differential equation Q is that 
the dynamics along (£1, £2) nas a simple diffusive behavior, since, by a straightforward 
Ito calculus, for any test function ip : T 2 — > R, 

9 t E(^((a,6)^)) =r 1 E(A^((ei,6)^)) (6) 

which is a weak form for the heat equation on the marginal law along (£1,^2) of 
the density of Y t . Roughly speaking, the energy landscape has been flattened in the 
(£1, ^-direction. 

The aim of this work is to explore a generalization of the ABF method — originally 
devised to compute a free-energy difference, which is a quantity of paramount impor- 
tance in statistical mechanics [14, 39_ — focusing primarily on its adaptive importance 
sampling feature and with the objective of obtaining a diffusion along some chosen 
directions. More specifically, we have in mind the case of m reaction coordinates with 
m > 4, for which the standard ABF approach cannot be used, because it would re- 
quire that the biasing forces, namely m functions of m variables, be approximated 
by a Monte Carlo procedure which is admittedly computationally prohibitive as m 
increases. Moreover, the fact that the biased dynamics along the reaction coordinates 
is a simple diffusion in the whole torus T m seems somewhat inappropriate, given that 
exploration of such a space may become extraordinarily long for large m. 

The approach proposed herein consists in considering the dynamics: 



dX t = - V(V - 2 A t &) {Xt) dt + ^2p^dB u 



a = l 



dA a 

for a = 1,2, ~^{x a ) = E{d Xa V{X t ) | £ a (X t ) = x a ). 



(7) 



The interest of this dynamics is that only two one-dimensional functions have to 
be approximated, ft is, therefore, expected that the Monte Carlo approximation of 
the biasing functions will be faster. One can check that this dynamics retains some 
essential features of the ABF dynamics p| , namely the fact that it leads to a simple 



diffusive behavior in each direction £ Q (see Section 2.1 below): for any test function 
tp : T -> R, and for a G {1,2}, 

d t E(p(UXt))) = ft' 1 E&"(UXt))). 

As a consequence, the marginal laws along £1 and along £2 of the equilibrium measure 
of the dynamics are uniform laws over T. ft ought to be noted, however, that ([6| does 
not hold in general in such a situation, and that the marginal law of the equilibrium 
measure along (^1,^2) is not in general a uniform law over T 2 . This is presented in 
detail in Section 12.11 

A motivation for considering dynamics ^ is that in the decoupled case, where 

V(x\, X 2 , X3.,. n ) = V(xi) + V(x 2 ,X3.., n ) or 
V(xi,X 2 ,X 3 ... n ) = V(x 2 ) + V(xi,X 3 ... n ), 

then (|7| is equivalent to Q. The key idea here is that, if the two reaction coordinates, 
£1 and £2 , are " n °t too strongly coupled", then ^ should be as effective as Q, at a 
far reduced cost. 

In Section [2] we propose a mathematical analysis of the long-time convergence 
of Q, which quantifies the underlying decoupling assumption. Section [3] is devoted 
to a discussion of some generalization of the idea of the present work, in particular 
to the case of m non- linear reaction coordinates with m > 2. Finally, in Section [4] 
we report two numerical illustrations on non-trivial test cases, which illuminate the 
interest and the limitation of the approach. 
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2 A convergence result 



Let us introduce the Fokker-Planck equation associated to Q. Let us further refer 
4>(t,x) to as the density of the distribution of X t . This function satisfies the partial 
differential equation: 

{ W = div (VV^ + /TV) - d Xl ({A\)'{xM) - d X2 ((A 2 t )'(x 2 )ij), 
d Xl V(x)tp(t, x) dx 2 dx 3 ... n 



K)>i) 



m(x 2 ) 



ip{t,x) dx 2 dx 3 ... n 
d X2 V(x)ip(t, x) dx x dx 3 ... n 
ip(t, x) dxi dx 3 ^ n 



(9) 



where (Af)' corresponds in what follows to the derivative of the one-dimensional 
function 

2.1 Diffusive behavior 

Let us introduce the marginal laws along x\ and x 2 of ip: 

ip Xl (t,xi)— J ip(t, x) dx 2 dx 3 ... n and ip X2 (i, x 2 ) = J ?/>(£, x) dx\ dx 3 ... n . (10) 
These marginal laws exhibit a simple diffusive behavior: 

Proposition 1 The probability distribution functions ip Xl and ip x ' 2 satisfy the heat 
equation: for a G {1,2}, 



P~ l d Xa , Xa V a =0 onT. 



(11) 



This property is easy to demonstrate by integrating the partial differential equation 
satisfied by ip in Q over the Xi, for i 6 {1, . . . n} \ {a}. 

As a simple consequence of (11 1, the marginal laws ip Xl and ip X2 converge to their 



equilibrium value It exponentially fast with rate 

r = 47r 2 , 

for example in the following relative entropy sense (see Definition [I] below) : 

(t, •) ln(^ (t, •)) < / f" (0, •) Hi> Xa (0, •)) exp(-2^- 1 rt). (12) 



2.2 Stationary state 

If A] and A\ reach a stationary state A^ and , it is standard that the stationary 
probability distribution function in ^ is: 

^(x) a eM-P(V(x) - AUxi) - A^x*))). 

Thus, proving the existence of a stationary state is tantamount to proving the existence 
of a couple (^4^ , A^ ) solution to (note that the functions A^ are defined up to an 
additive function): 



(^o)'(si) = 



d Xl V(x) exp(~P(V{x) - A^(x 2 ))) dx 2 dx 3 ... n 

exp(-j3(V(x) - Al (x 2 )))dx 2 dx 3 ... n 
d X2 V(x)exp(-/3(V(x) - A^xi))) dxx dx 3 .„ n 



(13) 



exp(-/3(V(x) - Alcixx))) dxx dx. 



3...TL 
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Let us set 



Finding a solution to ( 13 1 is then equivalent to find a couple (p 1 , p 2 ) solution to (note 
that the functions p a are defined up to a multiplicative constant): 



p\x 1 ) 



exp(-/3V(x)) 

P 2 {X2) 

exp(-/3V{x)) 



dx 2 dx 3 ... n , 
dxi dx 3 ... n . 



(14) 



Proposition 2 Let us assume that V is a continuous function on T n . Then, there 
exists a solution to (14), and thus there exists a stationary state (ip^, A^, A 2 ^) to ([9]). 



Proof : Let us build a sequence of continuous functions p\ : T — > as (pj = 1): 

l / \_rr\ f exp(~/3A(x 1 ,x 2 )) 
p n+1 (xi) - Z n+1 j dx 2 , 

/ rr~\ dx ^ 

where Z\ +1 is chosen such that J l//o^ +1 (xi) dx\ — 1 and A is the free energy ^ 
introduced above. It is clear that if (Pn) n >o converges in L°°(T) to p^, then p 1 = 

p^, and p 2 = [ — — — j- — ( Xl,X2 ^ dxi is a solution to (Tl4j) . We use the Arzela- 

J Poo{ X l) 

Ascoli theorem to show that (Pn)n>o is a compact sequence in the space of real- 
valued continuous functions over T, endowed with the L°°-norm. which concludes the 
existence proof. 

It first ought to be noted that since V is continuous, there exists positive reals a, b 
such that < a < exp(— f3A) < b on T 2 . We, hence, have, for all x\ G T, 

exp(-f3A(x 1 ,x 2 )) f exp(-pA(xi,x 2 )) , .a 
dx 2 > / = dx 2 > - . 



1 



exp(-pA(*i,X2)) dxi ^ J h f I " h 



Pn( X l) J Pnfcl) 

Likewise, 

exp(-f3A(x 1 ,x 2 )) , . b 
dx 2 < 



exp(-PA(xi,x 2 )) a 

1} \ dx i 

PnV x V 

>From this, one also obtains | < Z^ +1 < ^ which implies: Vn > 0, 

a 2 , b 2 

The equicontinuity property remains to be checked to conclude the proof with the 
Arzela-Ascoli theorem. This property is, however, straightforward to check by noting 
that 

b 2 f 

\pl{ x i) - pl{ x 'i)\ < Z2 / \ exp(~/3A{x 1 ,x 2 )) - exp(-l3A{x' 1 ,x 2 ))\dx 2 . 



This altogether proves the existence of a stationary state. Its uniqueness is a 
consequence of the convergence result stated in the next section, which holds under 
an additional weak-coupling assumption (see (|18[) below). 



A consequence of ( 14 1, which is also consistent with Qllh , is that the marginal 



probability density functions (see ( 1 1 ) of ip x along xi and x 2 are uniform: 



i>™( x i) = / ^oo(x) dx 2 dx 3 ... n = 1 and ip^(x 2 ) = / Voo(z) dx\ dx 3 ... n = 1. 
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Last, it is worth noting that by and large A^ and A 2 ^ are not the free energies 
associated to x± and x 2 and defined as — which should be compared with the definition 
of the bidimensional free energy |5| : 



A 1 (x 1 ) = -/J 1 ln (^J exp(-i3V(x l7 X2,x 3 ... n ))dx 2 dx 3 ... n ^j 
A 2 {x 2 ) = -/J -1 hi ^ exp(-f3V(x ll x 2 ,x 3 ... n ))dx 1 dx 3 .. 



(15) 



Actually, a special case for which = A 1 and A 2 ^ = A 2 (up to additive constants) 
is the decoupled case, namely if the two-dimensional free energy A (defined by (|5])) 
writes as a sum of a function of xi and a function of x 2 . which is equivalent, up to an 
additive constant, to: 

A(x 1 ,x 2 ) = A 1 (x 1 )+A 2 (x 2 ). (16) 
Under these premises, it is easy to check that p a = cxp(— (3A a ) is the unique solution 



to (14 1, up to a multiplicative constant. It ought to be noted that (K| implies (16) 



2.3 Convergence 

Let us now consider a stationary state (-000,^1^, ^4^,) to Q. The aim of this section 
is to prove the convergence of (|9| to this stationary state. 

Let us first introduce the conditional probability density functions: 

1po \x 1 \X2,X 3 ...n) = TxT/ \ = Voo(Xl 1 X 2 ,X 3 ... n ) 

Woo \ x l) 



I , ^ Woo\^l,^ 2 ,^ 3 ...n) , , x 

Voo\x 2 {Xl,X 3 ... n ) = ^ = 1poo{Xi,X 2 ,X3... n ). 



and 

(Xl 

r^{x2) 

In our particular case, this consists only in freezing one variable of ip^. 

To achieve this objective, we will need tools related to logarithmic Sobolev inequal- 
ities. We recall that (see [3l 141157]): 

Definition 1 A probability measure v is said to satisfy a logarithmic Sobolev inequal- 
ity with constant p > (in short: LSI(p)) if for all probability measures p such that 
p is absolutely continuous with respect to v (denoted p <i/ in the following), 

H{p\v) < ~I{p\v), 
2p 

where 

H{p\v) = J In dp 
is the relative entropy of p with respect to v and 

'dp 



l{p\v) = 



Vln , , 

av 



2 

dp. 



is the so-called Fisher information of p with respect to v . 

Since the measures 4> o\x 1 {x 2 ,x 3 ... n ) dx 2 dx 3 ... n and V'ooia^ (xi, Xs... n ) dx\dx 3 ... n are 
defined on a compact space (T™ _1 ) and since ip^ is smooth, it follows from standard 
arguments — e.g. Holley-Stroock criterion, that they satisfy logarithmic Sobolev 
inequalities. Let us introduce the associated positive constants p\ and p 2 : 

W~oo\x\ x 3 ... n ) dx 2 dx 3 ,,, n (resp. ipoo\x 2 (xi, x 3 _. n ) dx\dx 3 ,,,n) satisfies /i»t\ 
a LSI with constant p\ (resp. p 2 ) for all x\ e T (resp. x 2 € T). 
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We also need to introduce the coupling constants 

«1 = l|Va; :! ,a:3... n (9a!iV)||ioo( T «) and K 2 = ||V xl ,a !3 ... n (9 X3 V)||i a( T ») 

which are well defined since V is assumed to be smooth over T n . It is worth noting 
that «i = or ft 2 — is equivalent to (Jsl) , which implies the relation (16 1. This 
motivates the terminology of coupling constants. 

For the convergence result to hold, we need the coupling constants to be sufficiently 
small compared to the logarithmic Sobolev constants p\ and p?;. 

PlP2 > P 2 K1K2- (18) 

We are now in a position where we can state the main mathematical result of this 
contribution. 



Theorem 1 Let us assume (18 1. The probability density function tf>(t,-) then con- 
verges to ipao exponentially fast: for any e £ (0, A), 3C > 0, Vt > 



\ij)(t,x) — ipaa{x) \ dx < Cexp f — p 1 min ^(A — e),r^jt^j 



(19) 



where r = Air and 



\2 i 4h ik-2 



^ft+tt-yfo -«)' ■■ nn (20) 

is a positive constant. Furthermore, for any positive time to and e £ (0, A), 3C > 0, 
Vt > t , 

\( A t)' ~ i-Ko)'? <Cexp(- 2P~ X min ((A - e), rjtj , (21) 

where a € {1, 2}. 

The interpretation of this theorem is that, if the coupling constants k,\ and «2 are 
sufficiently small, the dynamics converges exponentially fast with a rate essentially 

limited by — ' — ^ (namely A when k\ = or k 2 = 0). This constant is expected 

to be larger than the logarithmic Sobolev constant of the original measure \i (which 
gives the rate of convergence of the original dynamics Q) if £i and £2 are well chosen 
— see related discussions in the work |38[ [35] . 

Proof : The proof is an adaptation of the proof for the long-time convergence 
of the ABF process, which can be found in [38 . It can be assumed without loss of 
generality that /3 = 1 up to the following change of variable: t = P~ 1 t, ip(i, x) — ip(t, x) 
and V(x) = PV(x). 

Let us first rewrite the partial differential equation satisfied by ip as: 



\ + d xi (((Air - {a])')4>) + o X2 (((Air ~ {Aiy)ii) . 



d t ijj — div I tp V 
Let us consider then the relative entropies 

E{t)=H{^ 00 )= f ^(t,.)ln(^(*,-)/^oo) 

and, for a € {1, 2}, 



The aim of the present proof is to show that E(t) converges to exponentially fast. It 
is already clear from (111 that Efj(t) converges to zero exponentially fast (see (12)), 
so that it is enough to consider 



E^{t) = E{t) - E%{t) = / ff(VvJV)IV>oc|*JV>* a (Ma)cfea, 
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where 



T"- 1 



is the relative entropy, with respect to ^oolxo.) of the conditional probability density 
functions 

, v i>(t,Xl,X2,X3...n) All* \ ' l P( t ' X l' X 2,X3...„) 

W\x 1 {t,x 2 ,x 3 ... n ) = ^ u and ip\ X2 (t,x 1 ,x 3 ... n ) = 



^(t, xl ) 



ip X2 {t,x 2 ) 



Let us focus on the case a = 1, albeit similar computations hold for a = 2. Let us 
compute 



dE x m _ dE dEjf 
dt dt dt 

Vln 



4>oc 



i> + Y, f ((^)'-(^o)')a*> 



Ipoc 



1> (22) 



\d Xl In^JlV 1 - 
It is rather easy to check the following identity: 



d x In 



dx^dx 3 ... n -d x Jn(ip x "), (23) 



where a = 1 (resp. 5=2) when a = 2 (resp. a = 1). Using (23 1 in (22 1, we obtain 



dt 



d Xl In 



V'oc 



V'oo 
2 



T VI" 



A) 

V'oo / 



?Ada;2 cfcc 3 ...„ ) — dx\ 



T" 



v 



d Xl m(V> xl )S xl ln( ^ )^ + / |9 xl ln(^)rV 



By virtue of the Cauchy-Schwarz inequality, the term on the second line is non-positive. 
Using again (23 1, we, hence, have 

dEl 



dt 



< 



9 x 2 ,x 3 ...„ In 



A 



if) - / d^lnir 1 )^ 1 ((AlY - (Al)') 



(24) 



We now need an estimate for |(.Aj)'— (A^)'!- For any coupling measure it € H(if)\ Xl (t, ■)> ^oo|a:i)j 
it holds: 

\(Aly(x 1 )-(AlY(x 1 )\ 

, X2j %3...n ) — d Xl V(xi,x' 2 ,x' 3 n )) n(dx 2 dx 3 .., n , dx' 2 dx 3 n ) 
\(x 2 ,x 3 ... n ) - (x' 2 ,x 3 n )\Tr(dx 2 dx 3 ... n , dx 2 dx 3 n ) 
\{x 2 ,x 3 ... n ) - (x 2 ,x 3 ^ n )\ir(dx 2 dx 3 ... n , dx' 2 dx 3 n ). 



/T™" 1 xT™ -1 

< \\V*i,x 3 ... n (dx 1 V)\\ L ° 



Taking the infimum over all ir G H({j, t , x , Moo, a;), we obtain 

\{A\)'{ Xl ) - «)'(z!)| < Kl W(^ Xl (t, 0,^1^ ) 
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where W stands for the (L 1 ) Wasserstein distance. We will now resort to the fact that 
if v is a probability measure satisfying a logarithmic Sobolev inequality with constant 
p, then we have the Talagrand inequality (see [HI HH]): For all probability measures p 
such that p <C 



W(jjl,v)< \l-HQi\v). 



Using (17 1 together with the Talagrand inequality, we, thus, obtain (the proof being 
evidently similar for a = 2): for a <E {1, 2}, 

\(A?)'(x a ) - (A^)'(x a )\ < ^-tfWvJV^cclxJ- (25) 

Using this estimate in (24 1, the constants pi and p2 introduced in (JTrj) , and Cauchy- 
Schwarz and Young inequalities, we get: 

2 



dEl < 1 
dt ~ 2 

1 

+ 2 



9.r 



hi 




ft'l 



-El 



P2 



Employing (111, it is standard to show that (see |38l Lemma 12] for example): 

I{r\t r )\^)<heM^rt) 
where I = /(V^ 1 (0, OIV'S ) an d> we recall, r = Air 2 . We finally find for any positive r/x- 



dEl 
dt 



■-El 



^/ cxp(-2r<). 



P2 " 2p\rii 

Utilizing a similar reasoning with a — 2, the following system of inequalities is ob- 
tained, wherein 771,772 are positive real numbers to be fixed: 



dE ln 

dt 

dE 2 m 
dt 

In the limit 771 = 772 = 



< -p^l-^El 

< - P 2{i-m)E 2 m 



P-2 



! -p i 



—^—I exp(-2rt), 
ZPiVi 



-EL 



Pi - 2p 2 i] 2 
0, we get the linear system: 

du 1 
~dt 
du 2 



Iq exp(— 2rt). 



< -piu 1 + -^-u 2 , 



Pi 



dt 



< 



-p 2 u 



Pi 



for which it can be be shown quite simply that, under the assumption (18 1, 
Vt>0, IKw 1 , it 2 )(*)|| < ||(M 1 ,u 2 )(0)||exp(-2Ai) 



where A is defined by (20 1. It is then easy to reach the result ( |19| , employing the 
Csiszar-Kullback inequality: 



(26) 



and the fact that H^ipoo) = E = E^ + E^. We refer the reader, for example, to 
the end of the proof of Theorem 1 in for a similar reasoning. 

Finally, the convergence results (21 1 on (A\ )' and (A 2 )' are easily obtained from (25 ) 
and the fact that (using (111) ip Xl and tp X2 are bounded from below by a positive con- 
stant for times larger than any arbitrary small positive time, see the beginning of 
Section 3.3.2 in [38 for more details. ^> 
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3 Discussion and generalizations 

3.1 The case of nonlinear reaction coordinates 

We now would like to discuss generalizations of the approach introduced above, in 
the case where (^1,^2) are n °t simply (xi,X2), which constitutes the vast majority of 
practical situations. 

Let us denote £1 : T> — > T, and £2 : 2? T, the chosen variables forming the 
multidimensional reaction coordinate. A natural generalization of ^ is the following: 

dX t = -V ( V - 2 A? o £ a ] (X t ) <ft + ^2^dB t , 

\ «=l / (27) 

for a = 1,2, ^-(*a) = E(/f (X t )|£ Q (X t ) = * a ), 
where /" stand for the so-called local mean forces associated to £ Q and defined by: 

/t_ l iv&i* " dlv l^aFJJ (28) 

and 

It ought to be noted that (27) reduces to ^ in the specific case of £ Q (x) = x a . 

If ^ denotes the probability density function of X t , then the marginal probability 
density functions are: 



which boils down to (10 1 in the case of £ a {x) = x a . 

The conditional measure 8^ar x \— Za (dx) obeys the following definition: For any test 
function tp : T> — > R, 



ip(x)dx= / / (p(x)5^^_ Za (dx)dz a , 
where E Q (z Q ) = {x G T>, £, a (x) = z a }. A corollary of the coarea formula is that: 

a -* £ ° = L M ( +div ^ ^w-'*'- (30) 



For detailed proofs and references, the reader is referred to [39, Lemma 3.10]. 



One can check the following property, reminiscent of (11 
Proposition 3 Let us assume that 

|V£i| = |V6| = I- (31) 

The probability distribution functions i/j^ 1 and tp^ 2 satisfy the heat equation: for a £ 
{1,2}, 

d t ^ - (3~ 1 d Za , Za ^ a = on T. (32) 



Proof : Let us prove (32 1 for a = 1, the proof being evidently similar for a — 2. For 
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any test function ip : T — > M: 
d t ip Sl <p = J d t ip<po 

= J div((V(V - A] o a - A 2 o £ 2 )^ + /3 _1 V^))<p o 6 

= - jf ((v(y - o^~a\o e 2 )^ + r ■ va ^ a 

= -/ v(v -A 2 t oe)-vz^v' + I / t W°£i 
= -r 1 / Ae^^'oei-r 1 / vv-va^°6 

Jv Jv 

= -p~ 1 [ d^V, 



T 



where we used (30 1 for the last equality. This is, indeed, a weak formulation of (32 1. 



Assumption (31) essentially amounts to assuming that £ Q is the signed distance to 
S a (0). It is straightforward to generalize the above result by changing the assump- 
tion plj ) to: |V£i| and |V^| are constant functions. It is also possible to generalize it 
to the case where |V£i| (resp. |V^2|) depends on x only through £i(x) (resp. £2(2;)) 
with slight modifications of the definition of the functions / t Q . Generalization of the 
convergence results of Theorem[l]to this setting is also possible, yet we will not pursue 
in this direction here. 

3.2 The multi-dimensional setting 

When more than two variables describing the multidimensional reaction coordinate 
are needed, the biasing procedure introduced above can be generalized as follows: 

(i) In the case of m collective variables £i,... , £ m , a natural generalization would 
consist in biasing the dynamics by a potential A\ o ^ + . . . + A™ o £ m , using 



a straightforward extension of the dynamics ( 27 1 . It is worth noting that the 
complexity is typically linear in m. whereas it is exponential in m for a standard 
ABF approach. For another biasing approach, in the context of multiple reaction 
coordinates, the reader is referred to |48| . 

(ii) In the case of m collective variables, one might want to keep certain collective 
variables coupled, e.g. £i and £2- A natural way to build a biased dynamics 
in such a case consists in considering a biasing potential A\' 2 o (£i,£2) + A^ o 
£ 3 + . . . + A™ o £ m , and using an adequate formula for updating A\ 2 based on 
the standard ABF approach Q (see |39l Section 5.1.1] for example for adequate 
formulae for general (£i,£2))- 

As mentioned above and roughly speaking, biasing the potential by the sum A\ o 
£1 + . • .+A™o£ m rather than by a m-dimensional adaptive potential A t o(£i, . . . , £ m ) — 
which would yield a perfect diffusion along the m-dimensional vector (£*, . . . , £ m )(X t ) 
- is tantamount to supposing some sort of decoupling on the collective variables 
£* , . . . , £ m . More precisely, if the free energy associated to (£1 , ... , £ m ) writes as a sum 
of functions A 1 o ^ + . . . + A rn o £ m , then the two algorithms (generalized- ABF and 
ABF) are equivalent. We, therefore, expect the method to be efficient if the collective 



11 



variables are loosely coupled, as would be the case for two dihedral angles distant 
from each other in a molecule. Conversely, should the two collective variables be 
more strongly coupled, it might be interesting to resolve this coupling in the adaptive 
potential, as discussed in item (ii) above. 

3.3 Further generalizations 



In practice, it may be difficult to implement the dynamics ( 27 1 , in particular on account 
of the computations of the analytical expressions for f" which may be cumbersome. 
A natural idea, which is, however, not supported by any mathematical reasoning, 
consists in simplifying the expressions for considering, for example: 

'-(W-'M^f))- 

or even 

f 1 = ^^P, (34) 



ive 



and similar formulae for f 2 . A very practical approach can be stated as follows: If 
A\ and A 2 happen to converge to some functions and A 2 ^, then it is always 
possible to fix this bias, and then to use unbiasing procedures akin to those described 



in Section 4.1 below — see (41 1, to obtain canonical averages 



4 Illustrations for free-energy calculations 

In this section, we illustrate the interest and the limitation of the approach in two test 
cases, using the Namd simulation package [47] . Propagation of motion is performed 
employing Langevin dynamics, in lieu of overdamped Langevin. Estimates of the 
biasing force rely upon trajectory, time averages rather than empirical averages over 



many replicas. Use is made of expression (33) to determine the local mean force. 
Overall, the method utilized here is similar to the biasing techniques ([7|-(p7| described 
above. 



4.1 Recovering the free energy 

To recover the free energy associated to some of the reaction coordinates chosen to 
bias the dynamics, one can simply use a grid of the values of the reaction coordinates, 
and classical formulas for the free energy or its derivative. This relies on the fact that 
from |7|-([27|, if ((Al)' , (A 2 )') reaches an equilibrium ((A 1 ^)' , (A 2 ^)'), then the law 
of Xt at equilibrium is proportional to exp(— j3(V — A^ o £t — A 2 ^ o £2))- 



Let us be more precise, considering (27 1. To get the free energy A(z\, 22) associated 



to the bidimensional reaction coordinate (^1,^2)) one can use the formula (compare 
with ([5]) which is (35) in the simple case (£1,62) = (^1,^2)): 



A(z 1 ,z 2 ) = -f3 1 ln / exp(-pV(x))5^ u £ 2 )( x )_( ZuZ2 )(dx). (35) 

It is indeed easy to verify that, at equilibrium, 

lim -/T 1 mE(<S e ((fc,6)(X t ) - *a))) + A^z,) + A^(z 2 ) = A( Zl ,z 2 ) (36) 

up to an additive constant. Here, S e denotes an approximation of identity converging 
to a Dirac mass at when e goes to 0. In practice, piecewise constant approximation 
is obtained over a grid of values accessible to the reaction coordinates. 

Another interesting formula to compute A(z%, z 2 ) derives from the formula: 

VA(z 1 ,z 2 )=E ll (F(X)\fa,&)(X) = (zi,z 2 )) (37) 
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where the notation E M means that X is distributed according to fi and F is a two- 
dimensional vector defined by: Va € {1,2}, 

F a = G a,7 V ^ ' VV - /T'div G- x 7 Ve 7 J , (38) 

7=1 \7=1 / 

where G~ ^ denotes the (a, 7)-component of the inverse of the matrix with components: 
Va, 7 e{l,2}, 

G an = V£ Q • V£ 7 . (39) 



It is easy to check that if the time marginal law of X t solution to (27 1 reaches an 
equilibrium, then 

E(F(Xt)\{Zu&){Xt) = C*i,*i)) = VA(«i,22). (40) 

At equilibrium, the law of X( is, indeed, proportional to exp(— /3(V— A^o^— A^ofa)), 
which only differs from fi by a multiplicative function of (^1,^2), which thus cancels 
out in the conditional average ( 40 ) . 

More generally, to estimate canonical averages, one may resort to the classical 
unbiasing procedure: 



L 



<p exp(-/3(^ ofr + A^o &)) Voo 
<pdl* = — r , (41) 

exp(-/3(A^ O^+^O &)) Voo 



•p 



where V'oo oc exp(— /?(V — o£j — A^ o£ 2 )) stands for the density of the equilibrium 
time marginal law of X t , once A\ and A\ have reached a stationary state. It is 
noteworthy that in practice, it is always possible to freeze A\ and A\ at a given fixed 
time to and apply the unbiasing procedure above with the bias A\ Q o £j + Af o o £ 2 
instead of A^ o^i+A^o^. 



4.2 Conformational equilibrium of the alanine dipeptide 

The first application of the method is a proof of concept making use of the prototyp- 
ical terminally blocked amino acid N-acetyl-N'-methylalanylamide (NANMA), often 
referred to as alanine "dipeptide" (51] . The molecular system consisted of NANMA 
immersed in a bath of 447 water molecules. Conformational sampling was performed 
over a period of 10 ns with the numerical scheme described here, in which the (j) and ip 
torsional angles of the backbone were handled as independent order parameters cov- 
ering each the [— 180;+180] range of the complete Ramachandran map [35] ■ In the 
present implementation of the method, the marginal biases, Al a ((j>) and A^ty), are 
periodical — i.e. the average of their derivative is expected to zero out. As a basis 
of comparison, a two-dimensional ABF calculation was conducted over a period of 
20 ns. For enhanced performances, the Ramachandran map was split into four in- 
dividual quadrants, corresponding to fully independent simulations. Furthermore, a 
standard molecular-dynamics simulation of equal length was performed, from which 
the 4> and ip dihedral angles were extracted to measure the exhaustiveness of the 
conformational sampling. 

The simulations were carried out using the Namd simulation package [47] in the 
isobaric-isothermal ensemble. The pressure and the temperature were fixed at 1 bar 
and 300 K, respectively, employing the Langevin piston algorithm |20| and softly 
damped Langevin dynamics. The molecular system was replicated in the three direc- 
tions of Cartesian space by means of periodic boundary conditions. The particle-mesh 
Ewald method [16J was employed to compute electrostatic interactions. The r-RESPA 
multiple time-step integrator [5B] was used with a time step of 2 fs and 4 fs for short- 
and long-range forces, respectively. Covalent bonds involving a hydrogen atom were 
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Figure 1: Conformational equilibrium of the alanine dipeptide (NANMA) in an aqueous 
solution. Distribution of the (<j>, ip) dihedral angles explored in the course of a 10-ns sim- 
ulation, wherein the torsional angles of the backbone were handled as independent order 
parameter (A). Reconstruction of the conformational free-energy map, employing expres- 
sion ( |36[ ) (B). Conformational free-energy map obtained from a 20-ns, two-dimensional ABF 
calculation (C). For comparison purposes, a similar map was generated from a 20-ns unbiased 
MD simulation (D). 



constrained to their equilibrium length. The short peptide and its environment were 
described by the all-atom Charmm force field [U] . 

In the past thirty years, the conformational equilibrium of NANMA has been 
investigated at different levels of detail, utilizing a variety of numerical schemes and 
potential energy functions |IIllin3[g|42l[50l[Il[^ 

I5T1 mil H51 H5] . Here, the investigation of NANMA is targeted at demonstrating the 
ability of the method to recover within reasonable sampling time the two-dimensional 
free-energy landscape that characterizes the conformational equilibrium of the peptide 
in an aqueous environment. As can be observed in Figure^ within 10 ns, the essential 
features of the ((f), ip) conformational space appear to have been explored. Plotting 
the number of samples accrued in the course of the simulation brings to light the 
anticipated coupling between the backbone torsional angles. Concomitant variation 
of the <p and ip angles is mirrored in the two parallel diagonals, which appear to 
be oversampled. After convergence of the marginal biases, A}^^) and A^^ifj), the 
two-dimensional free-energy landscape was reconstructed employing Equation ( 36 1 . 
Noteworthily, this map possesses the two expected pronounced minima corresponding 
to a right-handed a-helical conformation, often referred to as ckr, and to a /? strand, 
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together with ancillary local minima of higher free energy, associated to the so-called 
old and oll conformational states. 

Integration over the basins delineating the an and the /3 conformations yields a 
free-energy difference of about 0.2 kcal/mol, in favor of the former state, congruent 
with previous computer simulations [551 1501 IT51 HOI 1551 [T5] . Most importantly, the 
free-energy landscape inferred from the algorithm proposed herein is essentially iden- 
tical to that obtained from a 20-ns ABF calculation in ((f), tp) conformational space. 
A glimpse at Figure [T] is sufficient to conclude that not only the two maps appear to 
be almost interchangeable, but also, using the same bounds of integration over the 
otR and the /3 basins, the free-energy differences agree quantitatively. It can, however, 
be contended that the present toy-example may be somewhat exaggeratedly simple to 
be representative of more challenging instances, wherein the variables along which the 
local mean force is computed are more strongly coupled. As has been discussed above, 
this scenario would constitute a limiting case for the validity of the approach devel- 
oped in this contribution. Yet, as will be seen hereafter, estimators of the local mean 
force based on simulations handling unidimensional order parameters independently 
can still be used profitably to describe accurately rugged multidimensional free-energy 
landscapes wherein decoupling of the variables is not straightforward. Last, to illus- 
trate the role of importance sampling methods, the (0, ■0)-map regenerated from a 
20-ns unbiased molecular-dynamics simulation. It is apparent from Figure [T] that only 
the lowest free-energy states, i.e. otR and /3, have been visited, though sampling is by 
and large too parsimonious to allow an acceptably precise free-energy difference to be 
determined. 

4.3 Ion transport across a peptide nanotube 

In the second application of the method proposed herein, translocation of an halide 
ion through a chemically-tailored peptide nanotube is considered. Such engineered 
synthetic channels arise from the self-assembly through an intermolecular hydrogen- 
bond network [5] of cyclic peptides of alternated D- and L-chirality [231 12"5] , in which 
all the side chains are pointing outwards. Depending upon its amino-acid sequence, 
the resulting anti-parallel /3-sheet-like hollow tubular structure can further associate 
in the biological membrane with other channels to form nanopores, aggregate at the 
water-lipid interface prior to partitioning in the bilayer, disrupting in general the latter 
irreversibly, or simply span the membrane as an independent entity [22, 21 . Here, the 
peptide nanotube utilized consisted of eight stacked cyclic peptides containing alter- 
nated D-leucine and L-tryptophan residues and organized into «/do[LW]4 units [18_. 
At thermodynamic equilibrium, the cyclic peptides are 4.7 Aapart. The tubular struc- 
ture was immersed in a fully hydrated, thermalized palmitoylolcylphosphatidylcholinc 
(POPC) bilayer formed by 48 lipid units in equilibrium with 1,572 water molecules. 
The molecular assembly was replicated in the three directions of Cartesian space. The 
initial dimensions of the simulation cell were 36 x 41 x 79 A 3 . 

The free-energy landscape delineating permeation of a single chloride ion through 
the synthetic channel was determined along the longitudinal, £, and radial, p, direc- 
tions of the latter [57] . The surrogate, two-dimensional reaction coordinate (C,p) was 
constructed as a subset of cylindrical polar coordinates, namely the distance separat- 
ing the halide ion from the center of mass of the peptide nanotube projected onto its 
longitudinal axis, associated with the distance between the ion from this axis. Whereas 
the reaction pathway that connects the cytoplasm to the periplasm would span ap- 
proximately 40 A, the present investigation focuses on a 10-A segment, starting from 
the center of mass of the open-ended tubular structure. This restrained pathway was 
discretized in 0.1-A wide bins, in which force samples were accrued. The molecular- 
dynamics simulations were performed in the same conditions as described above for 
NANMA. 

Chemically -engineered peptide nanotubes possess the ability of conducting ions [52] . 
Assuming an appropriate amino-acid sequence, individual tubular structures can in- 
sert in the lipid bilayer, where they act as transmembrane channels [35]. That such 
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Figure 2: Permeation of a chloride ion through a peptide nanotube spanning a fully hy- 
drated lipid bilayer. £ denotes the longitudinal axis of the synthetic channel and p, the radial 



direction. Reconstruction employing expression ( 36 1 of the free-energy landscape of ion per- 



meation after 1 ns of sampling in which £ and p are treated as independent variables (A). 
Free-energy map obtained from a 1-ns two-dimensional ABF calculation (B). 10-ns ABF 
calculation, using as a starting point the two-dimensional gradients recovered from the 1-ns 
simulation in which £ and p are handled independently and employing expression ( |37[ ) (C). 
As a basis of comparison, (£, p) free-energy landscape inferred from a converged 30-ns two- 
dimensional ABF calculation (D). The lowest free-energy regions found at £ = —2.4 A and 
£ = —7.1 A correspond to an in-plane chelation of the halide ion, where the latter is located 
at the geometric center of the cyclic peptide. 



channels can be permeated by a small ions has been investigated at the theoretical 
level, employing a variety of methods O 031 HE] HI]- It has been suggested recently 
that whereas diffusion of a sodium ion through a synthetic channel formed of eight 
cyclo[LW]4 units and immersed in a fully hydrated POPC bilayer is essentially un- 
hampered [TS] , the same cannot be said for a chloride ion shuttled across the cavity of 
an identical peptide nanotube |27| . This result can be rationalized to a large extent 
by the lesser hydration number of the cation — viz. approximately 4-6, compared to 
that of the halide ion — viz. approximately 6-8 [5], which must undergo considerable 
dehydration to enter the open-ended tubular structure. In the midst of the latter, 
however, the free-energy landscape characterizing ion permeation is roughly similar 
for both species, in-plane coordination being enthalpically favored because it allows 
the ion to be better hydrated |18j . 
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As can be seen in Figure [2] this preference is reflected in the free-energy landscape 
obtained from the method proposed herein, handling £ and p as independent variables. 
Within 1 ns of sampling, the entire reaction pathway is explored, following what 
appears to be a minimum-action path. It is noteworthy that the chloride ion does 
not diffuse along a rectilinear path, collinear to the longitudinal axis of the synthetic 
channel, e.g. at £ = A, but rather avoids the large free-energy barrier of mid- 
plane coordination by grazing the wall of the tubular structure to interact with the 
amino groups of the cydo[LW]4 units. Repeated simulations, using different initial 
momenta, yield comparable landscapes, yet wherein the higher free-energy regions 
are essentially never visited. Interestingly enough, the free-energy minima found at 
about —2.4 and —7.1 A have virtually the same depth. It is remarkable that a limited 
simulation length of 1 ns would provide a consistent picture of the path followed by 
the chloride ion over the 10-A stretch of the peptide nanotube, when repeated two- 
dimensional ABF calculations of equal length only sample a small fraction of the (£, 
p) configurational space, as shown in Figure [2j 




10 -8 -6 -4 -2 -10 -8 -6 -4 -2 



<T(A) A C(A) 

Figure 3: Permeation of a chloride ion through a peptide nanotube spanning a fully hydrated 
lipid bilayer. £ denotes the longitudinal axis of the synthetic channel and p, the radial 
direction. Gradient of the free energy inferred from a 1-ns simulation wherein £ and p 
are handled independently, employing expression ( |37| (A). For comparison purposes, two- 
dimensional gradient of the free energy obtained from a 30-ns ABF calculation in (£, p) 
configurational space (B). 



Unfortunately, while the ABF calculation progressively explores the free-energy 
landscape of ion permeation, reaching convergence within 30 ns, a simulation of equal 
length based on the numerical scheme introduced here only marginally improves the 
picture drawn from the aforementioned short 1-ns run. This glaring shortcoming of 
the method can be ascribed in large measure to antagonist effects of the biases in 
the £ and in the p directions, deteriorating mutually any progress made by the two 
variables. It does not mean, however, that the valuable, yet incomplete description of 
the reaction pathway cannot be utilized profitably to recover the correct free-energy 
landscape, possibly faster than a classical ABF calculation would. 

As has been discussed previously, Equation (37) allows the gradient of the free 
energy to be estimated on the basis of independent measures of the force acting along 
the two order parameters, £ and p. In Figure [3] the reference two-dimensional gradient 
of the free energy inferred from a 30-ns ABF calculation is compared to an approxi- 
mation thereof obtained after 1 ns of sampling. Although the two vector fields show 
significant discrepancies, they also retain common characteristic features, in particu- 
lar in the region of lower free energy, visited appropriately by the algorithm described 
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herein. In turn, the approximate gradient can be employed as a starting point for 
a separate ABF calculation in (£, p) configurational space, with the hope that the 
initial guess of a minimum-action pathway might boost exploration of the complete 
free-energy landscape. As highlighted in Figure [2] this appears to be, indeed, true — 
after 10 ns, the map inferred from the separate two-dimensional ABF run possesses a 
topology essentially identical to that of the reference, 30-ns simulation. However not 
interchangeable, the two free-energy landscapes agree quantitatively in the low free- 
energy regions — i.e < p < 2 A, albeit only qualitatively so in the higher free-energy 
regions corresponding to the wall of the open-ended tubular structure. 
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